General recipe to realize photonic-crystal surface-emitting lasers with 100-W-to-1-kW single-mode operation

Realization of one-chip, ultra-large-area, coherent semiconductor lasers has been one of the ultimate goals of laser physics and photonics for decades. Surface-emitting lasers with two-dimensional photonic crystal resonators, referred to as photonic-crystal surface-emitting lasers (PCSELs), are expected to show promise for this purpose. However, neither the general conditions nor the concrete photonic crystal structures to realize 100-W-to-1-kW-class single-mode operation in PCSELs have yet to be clarified. Here, we analytically derive the general conditions for ultra-large-area (3~10 mm) single-mode operation in PCSELs. By considering not only the Hermitian but also the non-Hermitian optical couplings inside PCSELs, we mathematically derive the complex eigenfrequencies of the four photonic bands around the Γ point as well as the radiation constant difference between the fundamental and higher-order modes in a finite-size device. We then reveal concrete photonic crystal structures which allow the control of both Hermitian and non-Hermitian coupling coefficients to achieve 100-W-to-1-kW-class single-mode lasing.

It should be noted that the minus signs of the coupling coefficients in Eqs. (S2), (S3), (S4) ) and (4) when we re-define the coupling coefficients as follows; In Eq. (S14), θpc represents the phase of the non-Hermitian 180°-coupling (for example, Sx ⇒ radiative waves ⇒ Rx). θpc is determined by the position at which non-Hermitian couplings occur relative to the origin of the coordinate system (x=0, fixed at the center of the unit cell). For example, in a simple single-lattice photonic crystal whose circular air hole is placed at the center of the unit cell (x=0) (Fig. S1a), the position at which non-Hermitian couplings occur via a radiative wave (xr) is exactly located at x=0. In this case, the phase change in the ±180°-coupling θpc is equal to π because of our presupposition that the direction of the electric field vector is opposite for Rx and Sx. On the other hand, in a double-lattice photonic crystal shown in Fig. S1b, the situation is different: when the centroids of the air holes are equidistant from the center of the unit cell (x=0), the position at which the non-Hermitian couplings occur (xr) deviates from x=0. This is because the position x=xr is mainly determined by the center of gravity of the refractive index distribution, which is asymmetric with respect to the origin (x=0) in the double-lattice photonic crystal. Accordingly, the phase change of the ±180°-coupling θpc becomes 2βxr+π as shown in Fig. S1b, where β=2π/a. In our numerical analysis, θpc is almost constant, at 0.92π, because the variance of the hole parameters (d and 2x) are no larger than a few nanometers.
Since the value of θpc is dependent on the position of the air holes, its value changes following arbitrary global translations of the air holes as shown in Fig. S1c  Phase change of coupling coefficients following a global translation of the air holes along y=x.

Derivation of eigenfrequencies
From Eqs.
(2)-(5), we obtain pc pc pc pc Hermitian non-Hermitian non-Gamma Here, we assume that the deviation from the Γ-point (|kx|, |ky|) is much smaller than the reciprocal wavenumber β0. Consider the following matrix for basis transformation: 1 0 1 0 Block diagonalization of Eq. (S15) can be performed as follows; 1D pc pc pc pc 11 2D 1D 2D In Eq. (S17), the upper (lower) 2×2 block corresponds to modes A and C (B and D). The complex eigenfrequencies at the Γ point [Eqs. (9) and (10)] can be directly obtained from the first term of the right hand side of Eq. (S17). The eigenfrequencies of the modes in the Γ-M direction (kx=ky) can be also analytically solved, which gives the result shown in Eq. (16). When we consider the modes in the Γ-M' direction (kx=−ky), the coupled-wave matrix in Eq. (S17) cannot be diagonalized in a simple form. It should be noted, however, that the complex eigenfrequencies in the Γ-M' direction and those in the Γ-M direction are similar in the most cases owing to the following relationship: where I denotes a 4×4 identity matrix. According to this equation

Radiation constants of band-edge modes in a C2-asymmetric photonic crystal
As explained in the main text, when we consider a lattice-point design that completely breaks C2 symmetry, all of the band-edge modes in Eqs. (9) and (10) In one exemplary simulation, we change the complex coupling coefficient and fix the other parameters as follows: μ=90 cm -1 , 2D

Physical origin of radiation from photonic crystals and the role of the imaginary part I of the effective Hermitian coupling coefficient
To understand the physical origin of radiation from photonic crystals, we first review the role of the non-Hermitian coupling coefficient (iμ) in the radiation process.
Physically, iμ represents the vertical emission loss that accompanies each of the four fundamental waves as it couples to itself via radiative waves. This coupling coefficient iμ induces the radiation constant of 2μ for each fundamental wave (Rx, Ry, Sx, Sy) when the mutual coupling between them is ignored. For example, the radiation constant of each α A (cm -1 ) 0 180 Re (κ 1D + κ 2D-) (cm -1 ) Re (κ 1D + κ 2D-) (cm -1 ) On the other hand, at the Γ point of each band, the mutual couplings among the four fundamental waves determine the overall behavior of radiation from the four fundamental waves. As we have explained in Fig. 2 From this equation, one can easily understand that the coupling strengths from Sx + Sy to Rx + Ry and vice versa are different when the imaginary part of the phaseinvariant effective Hermitian coupling coefficients (I) takes a non-zero value. This asymmetry in the coupling strengths gives rise to the imperfect overall cancellation of radiation for mode A. When the condition of |I| << |R+iμ| is satisfied, the coupling strengths from Sx + Sy to Rx + Ry and vice versa are almost equal, which leads to a relatively small radiation constant of mode A (<20~30cm -1 ). This enables lasing with a low-to-moderate threshold current, which is desirable in a practical laser. To

Mode spectrum of a finite-sized PCSEL
In this section, we show an exemplary calculation of the mode frequencies and threshold gains (i.e., mode spectrum) of a finite-sized double-lattice PCSEL, and we demonstrate that the threshold gain margin between the fundamental mode and the higherorder modes originating from band-edge A corresponds to the "global" threshold gain margin, which considers the modes of all band edges, for single-mode lasing in large-area devices. We also discuss the symmetry of the electric field distributions of each mode in the finite-sized PCSEL. Figure S4a shows the mode spectrum of the 3-mm-diameter double-lattice Position for non-Hermitian couplings x r Node of electric fields of mode A Cancelled Not cancelled Not cancelled 13 PCSEL, which was designed for the simulation in Fig. 6 in the main text (R=45 cm -1 , μ=44 cm -1 , I=32 cm -1 ). From this figure, it is obvious that the modes originating from band-edges B, C, and D have much higher total loss than those originating from bandedge A; most of this difference is attributed to the much higher radiation constants of the modes of band-edges B-D, as explained in Supplementary Section 3 (note that the contribution of in-plane loss to this difference is much smaller, owing to the large area of the device). On the other hand, the difference in total loss between the fundamental mode (A-1) and higher-order modes (A-2, A-3) of band-edge A in a 3-mm-diameter PCSEL is much smaller than that between the bands, as shown in the inset. Therefore, it can be concluded that the threshold gain margin between the fundamental mode and the higherorder modes of band-edge A is also the global threshold margin of the device. Figure S4b shows the electric-field distributions of the modes with the smallest to fifth-smallest total loss in the same device. Here, we consider a photonic crystal which has reflection symmetry along y=x (Γ-M direction), and we assume that the shape of the current injection region is circular, which also maintains the reflection symmetry along y=x. As a result, the electric-field distribution of all modes in the finite-sized PCSEL become either symmetric or anti-symmetric along y=x. Specifically, the higher-order modes that have the lowest and second-lowest total loss (A-2 and A-3) have two antinodes in the Γ-M and Γ-M' directions, so their radiation constants can be described by simply 14 considering Δk along the Γ-M and Γ-M' directions. Other higher-order modes (A-4, A-5) have more complex electric-field distributions, but their threshold gains are much larger.
Therefore, in the discussion of the single-mode stability of PCSELs with reflection symmetry along y=x, it is sufficient to consider the wavenumber dependence of the radiation constants in the Γ-M and Γ-M' directions, which are considered in the main text.

Comparison between analytical formulae and numerical simulation of threshold margins of finite-size PCSELs
To verify the analytical formula of the radiation constant difference [Eq.
The simulated results of the radiation constant difference Δαv, in-plane loss difference Δα//, and threshold margin Δαtotal (=Δαv+Δα//) between the fundamental mode and the first higher-order mode are shown in Fig. S5 with red, grey, and black lines, respectively. As is apparent from this result, the radiation constant difference Δαv is much larger than the in-plane loss difference Δα// for 3-mm-diameter PCSELs.

Fourier component of double-lattice photonic crystals
As explained in the main text, the most dominant Fourier component that determines Here, Ntr is the transparency carrier density, gmax is the maximum gain, (−g0) is the absorption coefficient when there are no carriers, and ε is the gain suppression factor considering spectral hole burning and carrier heating. The detailed parameters used for the simulation are summarized in Table S2.
Taking the imaginary part of both sides in Eq. (S28), we obtain It should be noted that the right hand side of Eq. (S29) is constant. Considering the carrierinduced small change in Eq. (S29), we obtain Eq. (21) in the main text is directly obtained from Eq. (S31) in the vicinity of the Γ point

PCSELs
To experimentally realize the 100-W-to-1-kW PCSELs discussed in the main text, it is important to control the current density distribution of the device as well as to dissipate a large amount of heat generated inside of the device. Concerning the former, it is difficult to inject current into the center of a large-area device using a conventional ring-window electrode; to realize a spatially uniform current distribution, a mesh-type electrode, as was investigated in Ref. S9 for example, must be employed instead. In the case of a 3-mmdiameter PCSEL, by employing the mesh-type electrode shown in Fig. S6a, we can realize the current density distribution shown in Fig. S6b, in which the current density at the center of the electrode is only 20% lower than that at the edge. This current distribution is almost equal to that assumed in our time-domain simulations. It should be noted that we have numerically confirmed that the impact of the mesh on the beam profile is small because the width of each mesh (15 µm) is two orders of magnitude smaller than the diameter of the device (3 mm). In addition, the employment of the mesh-type electrode also contributes to the reduction of the series resistance of the device, which is useful to decrease the heat generated inside the device, which is discussed in the next paragraph. A similarly designed mesh-type electrode can be also applied to 10-mmdiameter PCSELs.
As for thermal management, the amount of heat generated inside of the device is negligible in the case of short-pulse operation with pulse widths of less than 1 µs, which is used for LiDAR and micro-processing. For continuous-wave (CW) operation, however, it is necessary to mount the device to a cooling package for heat dissipation. In this case, since the maximum heat dissipation per unit area is fixed, the maximum CW emission power is nearly proportional to the device area. Based on a CW emission power of ~7W 23 achieved using an 800-µm-diameter PCSEL [S10], it is expected that CW emission powers of ~98 W and ~1090 W are achievable using 3-mm-diameter and 10-mm-diameter PCSELs, respectively. For example, we have calculated that, for a 3-mm-diameter PCSEL mounted to a water-cooling package at 20℃ as shown in Fig. S6c, we can suppress the maximum temperature of the device to below 60℃ even under a heat generation of 200 W (which corresponds to an injection current of ~130 A and a CW output power of ~100 W) as shown in Fig. S6d. It should be also noted that the non-uniform temperature distribution shown in Fig. S6d and the resultant band-edge frequency distribution may affect the lasing characteristics of the device, but the effects of these distributions can be compensated by spatially varying the lattice constant of the photonic-crystal structure [S11]. Therefore, the 100-W-to-1-kW-class PCSELs proposed in this paper are experimentally feasible.